Ссылка на практикум
- Построим и "грубо" оптимизируем с помощью MOPAC структуры нафталена и азулена:
In [8]:
%%bash
echo 'C1=CC=C2C=CC=C2C=C1 Azulene' > azu.smi
obgen azu.smi > azu.mol
babel -imol azu.mol -omop azu_opt.mop -xk "PM6"
MOPAC2009.exe azu_opt.mop
babel -imopout azu_opt.out -opdb azu_opt.pdb
In [1]:
%%bash
echo 'c1ccc2ccccc2c1 Naphthalene' > napht.smi
obgen napht.smi > napht.mol
babel -imol napht.mol -omop napht_opt.mop -xk "PM6"
MOPAC2009.exe napht_opt.mop
babel -imopout napht_opt.out -opdb napht_opt.pdb
In [1]:
from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys
# pymol launching
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
### Если вывод в графическое окно тормозит или не нужен, то:
##__main__.pymol_argv = [ 'pymol', '-cp' ]
import pymol
pymol.finish_launching()
from pymol import cmd
### Информацию об ошибках можно смотреть там где запускали ipython notebook
### Будем вставлять файлы изображений
from IPython.display import Image
def Img(name):
cmd.ray()
cmd.png(name)
In [2]:
cmd.reinitialize()
cmd.load('azu.mol')
cmd.rotate(axis='x', angle=70)
cmd.do('set valence, on')
Img('azu.png')
Image(filename='azu.png')
Out[2]:
Азулен в конформации шезлонга
In [2]:
cmd.reinitialize()
cmd.load('azu_opt.pdb')
cmd.rotate(axis='x', angle=75)
cmd.do('set valence, on')
Img('azu_mop.png')
Image(filename='azu_mop.png')
Out[2]:
После преобразований Mopac выглядит плоским
In [7]:
cmd.reinitialize()
cmd.load('napht.mol')
cmd.rotate(axis='x', angle=-70)
cmd.do('set valence, on')
Img('napht.png')
Image(filename='napht.png')
Out[7]:
In [9]:
cmd.reinitialize()
cmd.load('napht_opt.pdb')
cmd.rotate(axis='x', angle=-70)
cmd.do('set valence, on')
Img('napht_mop.png')
Image(filename='napht_mop.png')
Out[9]:
С нафталеном глобальных изменений не произошло.
- Проведём более точную оптимизацию геометрии молекул с помощью GAMESS используя популярный базис 6-31G:
In [9]:
%%bash
babel -imopout azu_opt.out -ogamin azu_opt.inp
sed -i 's/ $CONTRL COORD=CART UNITS=ANGS $END/ $CONTRL COORD=CART UNITS=ANGS SCFTYP=RHF RUNTYP=OPTIMIZE $END\n $BASIS GBASIS=N31 NGAUSS=6 $end\n $system mwords=2 $end/' azu_opt.inp
babel -imopout napht_opt.out -ogamin napht_opt.inp
sed -i 's/ $CONTRL COORD=CART UNITS=ANGS $END/ $CONTRL COORD=CART UNITS=ANGS SCFTYP=RHF RUNTYP=OPTIMIZE $END\n $BASIS GBASIS=N31 NGAUSS=6 $end\n $system mwords=2 $end/' napht_opt.inp
In [10]:
%%bash
gms azu_opt.inp 1 | tee azu_opt.log
In [4]:
%%bash
gms napht_opt.inp 1 | tee napht_opt.log
- На основе полученных координат рассчитаем энергию методом Хартри-Фока и используя теорию функционала плотности:
Хартри-Фок:
In [13]:
%%bash
babel -igamout azu_opt.log -ogamin azu_hart.inp
sed -i 's/ $CONTRL COORD=CART UNITS=ANGS $END/ $CONTRL COORD=CART UNITS=ANGS SCFTYP=RHF RUNTYP=ENERGY $END\n $BASIS GBASIS=N31 NGAUSS=6\n POLAR=POPN31 NDFUNC=1 $END\n $GUESS GUESS=HUCKEL $END\n $system mwords=2 $end/' azu_hart.inp
gms azu_hart.inp 1 | tee azu_hart.log
In [5]:
%%bash
babel -igamout napht_opt.log -ogamin napht_hart.inp
sed -i 's/ $CONTRL COORD=CART UNITS=ANGS $END/ $CONTRL COORD=CART UNITS=ANGS SCFTYP=RHF RUNTYP=ENERGY $END\n $BASIS GBASIS=N31 NGAUSS=6\n POLAR=POPN31 NDFUNC=1 $END\n $GUESS GUESS=HUCKEL $END\n $system mwords=2 $end/' napht_hart.inp
gms napht_hart.inp 1 | tee napht_hart.log
DFT
In [14]:
%%bash
rm -rf /home/students/y12/popov/gamess-scratch/
babel -igamout azu_opt.log -ogamin azu_dft.inp
sed -i 's/ $CONTRL COORD=CART UNITS=ANGS $END/ $CONTRL COORD=CART UNITS=ANGS dfttyp=b3lyp RUNTYP=ENERGY $END\n $BASIS GBASIS=N31 NGAUSS=6\n POLAR=POPN31 NDFUNC=1 $END\n $GUESS GUESS=HUCKEL $END\n $system mwords=2 $end/' azu_dft.inp
gms azu_dft.inp 1 | tee azu_dft.log
In [6]:
%%bash
rm -rf /home/students/y12/popov/gamess-scratch/
babel -igamout napht_opt.log -ogamin napht_dft.inp
sed -i 's/ $CONTRL COORD=CART UNITS=ANGS $END/ $CONTRL COORD=CART UNITS=ANGS dfttyp=b3lyp RUNTYP=ENERGY $END\n $BASIS GBASIS=N31 NGAUSS=6\n POLAR=POPN31 NDFUNC=1 $END\n $GUESS GUESS=HUCKEL $END\n $system mwords=2 $end/' napht_dft.inp
gms napht_dft.inp 1 | tee napht_dft.log
- Извлечём значения энергий в энергетических единицах (Хартри):
In [6]:
%%bash
awk '/TOTAL ENERGY =/{print $4}' *_dft.log
awk '/TOTAL ENERGY =/{print $4}' *_hart.log
Учитывая, что \(E_{Hartree} = 627.509 469\; (kcal/mol)\), получаем таблицу:
| Вещество | E Naphthalene | E Azulene | ΔE, Hartree | ΔE, kCal/mol |
|---|---|---|---|---|
| Хартри-Фок | -383.3549 | -383.2822 | 0.0727 | 45.61994 |
| DFT | -385.6402 | -385.5852 | 0.055 | 34.51302 |
Из экспериментов известно, что энергия изомеризации нафталена в азулен составляет 35.3±2.2 kCal/mol. Согласно полученной таблице метод DFT оказался точнее метода Хартри-Фока, если судить по разнице энергий нафталена и азулена. При этом результаты DFT согласуются с экспериментом, в то время как метод Хартри-Фока дал ошибку ~22%.
- К сожалению, на kodomo доступ к Gabedit был запрещён, поэтому пришлось установить его себе локально.
В Gabedit было открыто окно "Display Geometry/Orbitals/Density/Vibration", в котором был загружен файл napht_dft.log. Затем была построена поверхность двух орбиталей.
HOMO (highest occupied molecular orbital):
34 E=-0.210800
In [2]:
Image(filename='HOMO.png')
Out[2]:
In [16]:
Image(filename='LUMO.png')
Out[16]:
Тоже самое для азулена
HOMO:
34 E=-0.186200
In [3]:
Image(filename='azuHOMO.png')
Out[3]:
LUMO:
35 E=-0.063300
In [4]:
Image(filename='azuLUMO.png')
Out[4]:
Сделать какие-либо выводы из этих картинок мне сложно, разве что исследуемые орбитали расположены на азулене и нафталене схожим образом, т.к. номера орбиталей HOMO и LUMO между этими молекулами совпадают.